
#### Kurtosis simulation ####
library(moments)

#### GSS analysis ####

gss_eqwlth <- read.csv("gssEQWLTH.csv")

sink("gssEQWLTHdescriptives.txt")
mean(apply(gss_eqwlth[,-1],2,sum))
max(apply(gss_eqwlth[,-1],2,sum))
min(apply(gss_eqwlth[,-1],2,sum))
length(apply(gss_eqwlth[,-1],2,sum))
sink()

eqwlth <- NULL
year <- NULL
yrs <- as.numeric(substr(names(gss_eqwlth),2,5)[2:ncol(gss_eqwlth)])
for(i in 2:ncol(gss_eqwlth)){
  for(j in 1:7){
    obs.num <- gss_eqwlth[j,i]
    eqwlth <- c(eqwlth,rep(j,obs.num))
    year <- c(year,rep(yrs[i-1],obs.num))
  }
}


wealth_year <- data.frame(eqwlth,year)

observed.variances <- c(by(eqwlth,year,var))

nboot <- 1000

set.seed(1234)

boot.variances <- matrix(NA,nboot,length(observed.variances))

for(b in 1:nboot){
  boot.index <- sample(1:nrow(wealth_year),nrow(wealth_year),rep=T)
  boot.data <- wealth_year[boot.index,]
  boot.var.b <- by(boot.data$eqwlth,boot.data$year,var)
  boot.variances[b,] <- boot.var.b
  if(b/100==round(b/100)) print(b)
}

boot.ci <- apply(boot.variances,2,quantile,prob=c(0.025,.975))

pdf("gssEQWLTHvariance.pdf",height=4,width=6)
par(las=1)
plot(yrs, observed.variances,type="n",ylim=c(min(boot.ci[1,]),max(boot.ci[2,])),ylab="Variance",xlab="")
abline(h=seq(3.4,4.6,by=0.2),v=seq(1978,2016,by=2),lty=3,col="grey55")
segments(x0=yrs,x1=yrs,y0=boot.ci[1,],y1=boot.ci[2,],lwd=3,col="grey40")
points(yrs, observed.variances,pch=19)
title(xlab="Year",line=2.25)
dev.off()


# BP test for heteroskedasticity
sink("gssEQWLTHbptest.txt")
library(olsrr)
model <- lm(eqwlth~as.factor(year))
ols_bp_test(model,rhs=T)
sink()

# series kurtosis analysis
## overall item distribution
table(wealth_year$eqwlth)
## kurtosis of the differences in means
kurtosis(diff(by(wealth_year$eqwlth,wealth_year$year,mean)))
## maximum within-year kurtosis
max(by(wealth_year$eqwlth,wealth_year$year,kurtosis))
## mean within-year kurtosis
mean(by(wealth_year$eqwlth,wealth_year$year,kurtosis))




